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Preface 


The  itudy  fCfKMMd  hoclB  wai  auchorued  by  HeadqoKicnt.  U  S.  Anay 
Coffw  of  Eagiaien  (HQUSACE).  «  pen  of  che  HyAaulict  Probkn  Aice  of 
(he  Repeif .  Evelmdoa.  Maiwiwieiife.  nd  fUhabitiiaQoe  (REMR)  RcMserdi 
Profram.  The  work  wee  perfenned  tnder  Ctvtl  Works  Research  Work  Uak 
32657,  ’Model  for  Evkhieiaoa  aad  Meiieenenre  of  Hlfh-Vdociiy  Channels.* 
The  REMR  TcdMMcal  Monitor  wm  Mr  Dive  Wh^  (CSCW-EH) 

Mr  WiUtam  N  Rnahtu  (CERD-O  wae  ihe  REMR  Coordineior  m  the 
Ouecioraie  of  Reseerch  and  DevdopaMni.  HQUSACE  Mr.  lanaeB  E.  Crews 
(CECWO)  and  Or.  Tony  C  Let  (CGCW-EG)  aarwed  as  the  REMR  Overview 
Coennunee  The  asrsaa  Ptoeraui  ManMer  was  Mr  WiitiBHi  F.  McOeeae. 

U  S  Army  Eafiaeer  Waterways  EaperaaMi  Station  (WES).  Mr.  OlennA. 
Pkkermg.  Chief.  Hydranlk  Straemres  Division  (HSD).  Hydrmdics 
Laboratory  (ML).  WES.  was  dn  Probleni  Area  Landv  Dr.  R.  C 
Berfer.  Eatuariae  Dmuon  (ED).  HL.  and  Mr  Rirtiard  L  SaodocUI.  H^. 
were  Friacipa)  btveatifaiors 

That  work  was  condnoad  by  Or  Boper  and  Mr.  Stockarffl  nnder  ihe  tfrect 
Mipervieioo  of  Mr  Pkkarim  ^  fsnaral  aayarviaion  of 

Messrs  a  Hcrnaann.  Ir  .  Diracmr  HL;  Richard  A.  Saner.  Aaristaai 
Director.  ML.  WiOiani  H.  McAnaHy.  Ir..  Chief.  ED;  and  John  P.  Oeotfe. 
Chief.  Locks  Md  Coadmis  Brandt.  KSD.  ML  Mr  Mike  W  Ob.  HSD. 
prepared  the  UhsKraiioni 

Ai  Ihe  liBK  of  pnMkaiion  of  dds  repon.  Direcaor  of  WES  was 
Or  Roben  W  Whalin.  Coawawtw  was  COL  Brooe  K  Howwd.  EN 
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Conversion  Factors, 
Non>SI  to  SI  Units  of 
Measurement 


Noe'S!  ttotti  of  (BCMwcmcM  Mcd  im  Si»  n^oi%  cm  tie  ooevcfietf  id  SI  mim 
m  Mkfwt: 


1  Introduction 


Background 

Tbe  hy^riulic  pgrfonmiicc  of  » iMflhvdociiy  fHwwfl  depends  on  immwm- 
iag  a  tupocriiical  Aow  rcgia»  ovcf  ipnciltod  pevttOM  of  in  Imgiti 
Predictifif  the  p«*t**^  of  ihocke.  ae  <yHiH|ttir  f?— wsvee  end 

hydraulic  junpa.  Md  the  wpcrckvation  of  Che  wtm  mrtact  u 

channel  benda  arc  neccnary  lo  cvahaMc  and  meiniain  the  icqewl  '■**11 
height*  Typically  ewpirkal  agpurtone  or  phyikal  hydraulic  model*  have 
been  used  u>  make  theac  evahmiion* 

Physical  model*  were  wed  in  the  otigmnl  «edy  of  many  caianni  Aood> 
control  channeli.  but  urhaniiatinn  whhln  thetf  drainege  baaine  ha*  remhad  in 
disclairgc*  greater  daw  thoee  for  which  the  chmwei  was  origiInhDy  daugned- 
Obaiaclea.  such  a*  dahri*  or  bridge  piers,  amy  cauw  the  flow  lo  lump  to  a 
subcritical  Mat*,  thw  reaulling  in  flood  dwaage  Therefore,  m  awapnstee 
and  a  readily  svadaMc  BMaa*  for  cvahHiiag  ihaae  channeli  is  needed.  A 
oumerical  modri  ia  a  logicai  approach. 

A  manerical  flow  modri.  HIVEL20.  has  bane  dtteeloped  a*  a  tool  to 
cvahiaca  hi0-«Uociiy  rhatniel*  which  ant  nanMnadr  ooaciwe  Unad  channel* 
with  hydraolkatly  ateqp  slopw  The  succaashd  dasign  of  Wgh-eelocicy 
channeli  depends  on  acconne  predketona  of  the  flow  depth*.  HIVELIO  i*  a 
depth-amaged.  rwo-dhnenBioaal  (3<0)  flow  awdal  diul|pMid  ipaeflkaily  for 
flow  fluid*  that  contain  wperaritical  and  whcrlflcal  reghnaa  a*  weO  as  the 
iraneitione  between  the  reghnee  The  model  doee  aoi  inUade  CorioUs  or  srind 
•ffocts  and  — ti— irmawnt  is  nqf  eddteBMd  since  dtme  nhenoBMOS  ere  not 
applicable  to  concrete  tiacd  chennalB 


Purpose  Slid  Scopo 

The  purpoee  of  ihii  report  ie  to  dcacribe  dm  naaHtriril  flow  model. 
HIVELID.  aad  to  Uhanraic  typkhl  Ugh-eefodey  flow  field*  that  the  modd  is 
capiMc  o€  linBNiating.  Thwa  modd  eeriflcaiioai  conU*t  of  oowparlng  renali* 
compiiiad  whig  HIVEL20  whh  Idboraiory  dma.  Modd  thniiaiioa*  are  ateo 
disewsed.  h  h  haporw  that  HfVeLlO  wen  endetstand  whfch  flow 
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situatkms  are  suited  to  simuiatioa  using  HIVEL2D  and  wliich  are  not. 

A  user’s  manual  will  be  the  second  report  ptdilished  for  this  work  unit. 
The  user’s  manual  will  discuss  such  topics  as  grid  generatioo  (FastTABS), 
hardware  requireraems  (personal  computers  (PC’s)  and/(H  workaatiaos).  time 
uivolved  in  grid  genenuioo  and  simulation  runs,  and  piidance  relative  to  grid 
layout  and  parameter  chotce 
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1  Inirodwctaon 


2  Description  of  Modei 


Tbe  model  is  designed  lo  iwmtisir  flow  typicitl  in  h^vHocity  ditnneit. 
The  model  is  s  fimie  ekmem  dcscrigiioe  of  dm  2-D  shaUow-wiicr  cquttkws 
in  conservative  form.  Specdkally.  die  model  is  designed  lo  provide  the  water 
surface  in  and  around  boundwy  transtetons.  bridge  piers,  eonltawices.  bends, 
and  other  gcoaccric  feamres  common  m  btgh-vdodiy  channels.  HIVEL2D  is 
^icabk  to  tub-  and  supercritical  flow  and  the  twnsatiews  bmween  dwse 
regimes 

The  sviicm  of  nontinem  stsaMioos  arc  solved  using  dm  Newton-Rsphson 
iterattve  method  The  Hewioii-IUplHOtt  medmd  was  sclscted  for  ihii 
bccauM  the  nonlinear  icrms  in  h^vdoeny  channel  flow  flelds  are  <puie 
sifmfkani 

Stresses  arc  modeled  ma^  dm  Mamung  s  formutemn  for  boundary  drag 
and  dm  Bousstnest) rclaiion for  Reynolds  totmei  Eddy  visoosioesme 
approaunaicd  using  w  empirical  rdatton  basod  on  Meming's  oocffkmnt  and 
local  flow  vsnabks 


Assumptions 

HrvEUO  solves  dm  depth-averaged  unsieady  shnDow  waier  nqpmti^ 
unplteiily  usM«  flnim  daneMs  The  governing  e^iaiiona  are  presented  m 
detail  in  Appendix  A 

An  unDotiani  (ssid  scanted)  aseueBpeioo  sadc  in  dm  dnivadon  of  dm  gov¬ 
erning  tqmkm  ia  dmi  vertical  acceteations  are  oegtigihle  when  compared  to 
the  horraomd  accdcratkma  and  aoedeteion  tee  »  gnivtty  Thd  assumes  thss 
the  pressure  (hstdbution  is  hydrustaiic  Tbe  hydrosiaiic  premure  maurapdon 
is  gcimraDy  appIkaMc  to  h^i-vdociry  ebawmh  However,  cm  mua  be 
taken  in  applying  dm  modd  to  flow  sicuteom  where  dm  wavUengfo  is  not 
long  reteive  to  the  depdi.  All  hydroeiacic  modeb  urcrsMinMc  the  speed  of 
short  wavdengdis 

To  undBimnd  why  short  wsvdengdts  trrrd  more  ilowfy  dwn  kmger 
wavetci«ihi.  consider  the  veiticd  aomenwm  erpanioo 


oweMr  i 


dW 

~S 


(1) 


1  dP 

.  m 

P  dz 


wbere 

p  »  fluid  denciiy 
F  «  preuure 

z  «  vcftkal  Cartesttn  coocdinaie 
w  »  the  z-diractkM  oonyonent  oi  veloctty 
t  •  tine 

g  »  eccekmioa  due  to  gravity 

Near  the  crest  of  s  wave,  or  anywhere  that  the  water  surface  is  convex,  the 
fluid  ekmot  rkhng  near  the  sitffncc  has  a  negative  vertical  acceleration,  i.e.. 

<  0.  Therefore,  the  vertical  acockrarion  acts  opiwaite  to  gravity  and 
pressure  is  reduced.  Convendy.  near  the  trough  where  the  water  surkoe  it 
concave,  the  oear-eurfKe  fluid  dement  has  a  positive  accekration.  and  to  acts 
to  enhance  the  pressure  The  net  effect  it  that  the  surface  curvature  reduces 
the  pressure  gradkott  and  the  enve  will  travd  mote  slowly.  The  shorter  the 
wavekngih.  the  slower  it  will  tend  w  propagate.  For  long  wavekngths  this 
unpact  ts  negl^ibk. 

Another  consequence  of  the  hydrostatic  assumption  a  that  the  ahallow- 
w«ar  cquattona  can  provide  no  kforakion  about  the  length  of  a  jump,  i.e., 
disunce  frotn  the  leading  10  the  trailing  edge  of  a  hydraulic  jump.  The 
shallow'Waier  enuttioni  no  vertical  velocilv  or  aocekratioo;  therefore, 
any  energy  that  ihould  he  captured  in  vertical  motion  is  lost.  The  shallow- 
wamr  equationa  mutt  stmukie  the  jump  aa  a  diaconrinuiiy  mid  all  vertical 
energy  is  wnmadiatdy  dkaipaiad  Thk  is  unIke  the  real  tytiem  in  educh  an 
umhikr  pnnp  may  farm  that  diasipaiei  tltis  vertical  motion  over  a  long 
distamx 

Another  signifleam  asaunguion  k  dmi  the  bed  slope  «  geomrtrkalty  mild 
sothatsinfotHi*#.  where  #  k  the  channel  dope  Tltia  geoaccrically 
ouid  slope  aaswnprion  k  dktingtikhed  from  a  bydkntically  mtid  slope  that 
conveys  flow  m  subcriticai  (i  c  .  the  timilow  wamr  equationa  assume  that  the 
slope  k  feamcerkafly  ntiM.  although  it  may  he  hydrautkany  amp).  Thk  k  a 
reaaonaMe  aaaumption  in  moat  high-velociiy  channd  applkitioni  where  the 
siope  k  ten  dun  0  02  However,  sehen  appikkion  k  mnde  to  a  long  channri 
reach  having  a  kvorabte  skpc  in  excem  of  say  O  OS.  dw  mild  dope 
Msumption  snM  tend  to  ovrrwfimBit  die  flow  qpeed  and  underterimaap  the 
Row  <kpth. 
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Discretization 


Discrete  values  of  the  unknown  variables  are  solved  using  a  Petrov- 
Galerkin  finite  element  formulation  of  the  governing  equations.  Details  of  the 
Finite  element  formulation  are  presented  in  Appendix  B.  The  model  is  quite 
stable  as  a  result  of  the  particular  Petrov-Galerkin  lest  function  employed, 
which  is  weighted  upstream  along  characteristic  lines.  The  finite  element 
model  reproduces  triangular  and  quadrilateral  bilinear  elements  and  can  be 
easily  modified  to  capture  complex  geometries. 

HIVEL2D  solves  the  finite  element  equations  using  the  weak  formulation. 
The  weak  formulation  facilitates  the  ^lecification  of  natural  boundary  condi¬ 
tions  and  allows  accurate  shock  capturing,  which  is  essential  in  modeling 
supercritical  flow. 

Temporal  derivatives  are  solved  implicitly  using  a  fmite  difference  discreti¬ 
zation.  The  temporal  derivatives  can  be  approximated  using  either  a  first-order 
(a  =  1.0)  or  second-order  (a  =  1.5)  backward  difference.  A  First-order  differ¬ 
ence  is  used  for  the  spin-up  to  a  steady  flow  solution,  whereas  a  second-order 
difference  is  more  appropriate  for  unsteady  flow  simulation. 


Boundary  Conditions 

Inflow  and  outflow  boundary  conditions  can  be  specified  as  supercritical  or 
subcritkal.  Supercritical  inflow  boundaries  require  the  qpeciFication  of  the 
three  degrees  of  freedom  (i.e.,  the  x-  and  y-components  of  unit  discharge  and 
the  depth  or  the  x-  and  y-components  of  velocity  and  the  depth).  Subcritical 
inflow  boundaries  i  |uire  that  only  two  ^leciFic  degrees  of  freedom  be  given 
(i.e.,  the  x-  and  y-components  of  unit  discharge  or  the  x-  and  y-components  of 
velocity).  Bocndary  conditions  at  supercritical  outflow  boundaries  should  not 
be  specified.  Subcritical  outflow  boundaries  require  the  specification  of 
lailwater  elevation  (i.e..  depth  plus  bed  elevation). 

Sidewall  boundary  conditions  are  enforced  using  the  weak  statement  The 
boundary  conditions  enforced  are  that  the  mass  and  momentum  flux  through 
the  sidewalls  are  zero.  Sidewalls  also  enforce  a  partial  slip  condition. 


OactwZ 


3  Test  Cases 


Numerous  test  esses  are  presented  to  illustrste  the  model’s  ability  to  simu¬ 
late  higb-velocity  cbanoei  flow  fields.  More  iinfx>itantly.  results  of  test  cases 
point  out  flow  conditioas  that  are  not  accurately  modeled  by  HIVEL2D.  The 
model  user  must  understand  when  the  model’s  governing  equations  are  ade¬ 
quate  for  the  patticular  problem  to  be  solved  and  when  they  are  not.  It  is 
(kmonstrated  that  the  numerical  scheme  accurately  solves  ^  depth-averaged 
2-D  equations.  The  depth-averaged  equations  appropritudy  describe  most 
higb-velocity  chatmd  flow  fields. 

Several  fetfures  conanonly  found  in  high-veloctty  channels  are  included  in 
the  test  cases:  simulatioos  of  channel  contractions,  expansions,  confluences, 
bends,  hydraulic  jumps,  bridge  piers,  bora,  and  toll  waves. 

Input  parameters  for  each  simulation  are  provided  in  tabular  form:  the 
temporal  difference  coefficient  a;  the  Petrov-Galerldn  weight  coefficieni  $ 

(for  smooth  flow)  and  dj  (for  flow  near  shocks);  Manning’s  coeffici  mt  n;  the 
eddy  viscosity  p  (for  smooth  flow)  and  v,  (for  flow  near  shocks);  and  the 
acceleration  due  to  gravity  g.  Detailed  deKtiptioos  of  the  temporal  difference 
and  Petrov-Oalerfcin  weight  coefficients  are  provided  in  Appendix  B.  Early 
tests  were  conducted  using  a  verskm  of  HIVEL2D  that  used  a  connant  t»er- 
specified  eddy  viscos^  for  the  entire  flow  field.  The  current  verskm  of 
HIVEL2D  varies  the  viscosky  baaed  on  local  flow  variables.  Details  of  how 
viscosities  are  determined  are  presented  in  Appendix  A. 


Channal  Contraction 


Initially.  HIVEL2D  was  run  to  simulate  supercritical  ffow  m  a  channd 
ooikractkm.  A  labwatmy  invesUgatkm  of  supercrkical  flow  in  chamid 
contractkms  was  reported  by  Ippen  and  Dawson  (1951)  m  which  depth  con¬ 
tours  were  presented.  One  particular  teat  doctaoeittd  flow  condkionB  in  a 
straight-wall  contnikioo  with  an  approach  Froude  number  of  4.0.  The 
straight-wall  contractioa  was  a  2-ft-  (0.61-m-)  wide'  chamd.  transitioning  to 


'  A  ttbte  of  fiKlon  for 
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a  I'ft-  (0.30-11)-)  wide  chaiuiel  at  a  convergence  angle  of  6  deg.  The  grid 
used  for  the  simulation  is  provided  in  Figure  1.  The  grid  consisted  of  1661 
nodes  and  ISOO  elements.  The  long  ap(m)ach  length  of  20  ft  was  used  to 
ensure  uniform  flow  conditions  at  the  ctmtraction.  Figure  2  shows  contours 
of  flow  depth  observed  in  the  laboratory  and  those  comiMited  using  HIVEL2D 
for  a  discharge  of  1 .44  cfs  (0.041  cms)  and  an  upstream  Froude  number  of 
4.0.  li^t  parameters  used  in  the  numerical  simulation  are  provided  in 
Table  1.  A  water-surface  mesh  generated  using  these  HIVEL2D  results  is 
shown  in  Figure  3. 

The  comraction  is  a  geometrically  simple  case,  but  the  results  demonstrate 
the  model's  ability  to  capture  the  oblique  standing  waves  resulting  from 
changes  in  the  wall  boundaries.  Oblique  standing  waves  are  created  as  the 
supercritical  velocity  encouttfers  the  converging  sidewalls.  These  standing 
waves  are  reflected  from  opposing  sidewalb  for  a  distance  downstream. 
HIVEL2D  captures  the  standing  waves  and  the  depth  increase  as  the  waves 
intersect  at  the  channel  center  line. 


Tahto  1 

Input  Pairamatars  for  tho  Channof  Contraction 

CowShtow 

VUm 

a 

to 

0.1. 0.5 

n 

0.0)1 

0.005.0.006  ft*/*«c  (4.86  *  (O'* 

0 

32.208  tXJ9»d  (9.817  m/MC*l 

Hydraulic  Jump 


Tests  were  conducted  to  compare  the  results  from  H1VEL2D  with  some  of 
those  obtained  in  the  physical  model  study  of  Rio  Puerto  Nuevo  Flood-Control 
Channel  (Stockstill  and  Leech  1990)  conduaed  fm  the  U  S.  Army  Engineer 
District.  Jacksonville.  Three  models  were  constructed  and  tested  for  the 
Puerto  Nuevo  study.  These  models  were  seleaed  for  comparisons  with 
HrVEL2D  calcultfiom  simply  because  the  channels'  geometric  and  hydraulic 
details  were  at  hand  and  several  features  common  to  high-velocity  channels 
were  included.  Speciflcally.  the  Puerto  Nuevo  study  channels  contained 
supercritical  confluences,  horizontal  curves,  expansions,  and  a  tailwater- 
imposed  hydr»)lic  jump.  To  avoid  similitude  questions,  HIVEL2D  was  run  at 
laboratory  scale.  his  resulted  in  direct  comparisons  of  calculations  to 
observed  quantities. 

The  Margariu  Channel  was  one  of  the  three  models  tested  for  the  project. 
The  Margariu  Channel  consi^ed  of  two  reverse  curves  separated  by  a  short 
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tangent,  a  channel  expansion,  and  Utiwater  that  produced  a  hydraulic  jump 
upstream  of  the  width  transition.  Details  of  the  modeled  rea^  of  the 
Margarita  Channel  are  preseitted  in  Figure  4.  Ir^nit  parameters  used  in 
UKxleling  the  Margarita  Channel  are  provided  in  Table  2,  and  the  fmite  ele¬ 
ment  grid  is  shown  in  Figure  S.  Typically  when  steady-state  solutions  are 
sought,  a  first-order  temporal  derivative  is  appropriate  and  an  a  =  1.0  is 
used.  Because  portions  of  the  Margariu  Channel  exhibited  unsteady  flow 
feanires,  an  a  »  l.S  was  used  resulting  in  a  second-order  temporal  derivative. 
The  finite  element  grid  consisted  of  796  nodes  and  6S1  elements.  Figure  6 
compares  water-surface  profiles  along  the  channel  walls  for  a  discharge  of 
1 .2  cfs  (0.034  cms)  with  a  tailwater  depth  at  the  downstream  end  of  the  model 
of  0.601  ft  (0.183  m). 


Tabl«2 

Input  Paramatart  for  tha  Margwfta  Charmal 

ConMon 

VakM 

a 

l.S 

AA 

0.2,05 

n 

0.010 

V.¥, 

0.008.0.026  ft^/aac  17.43  *  10^2.32  x  lO^mVsac) 

8 
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HIVEL2D  not  only  captured  the  hydraulic  junq>  in  the  channel,  but  also 
reproduced  the  asymmetric  flow  patterns  downstream  of  the  width  transition 
that  were  observed  in  the  physical  model.  The  short  length  of  the  transition 
( I  transverse  on  4  longitudinal)  in  conjunction  with  the  asymmetric  flow 
distribution  produced  by  the  upstream  bend  resulted  in  a  large  eddy  in  the 
downstream  channel  thtf  produced  flow  concentrations  along  the  left  wall. 
Depth-averaged  velocities  at  three  stations  are  presented  in  Figures  7-9. 

Higher  velocities  (1.9  fps  (0.58  n^M))  existed  ^ong  the  left  wall  whereas  the 
flow  along  the  right  wall  was  essemially  sugnant  as  shown  by  the  dark  area  in 
Figure  10.  The  numerical  model  results  shown  in  Figure  11  show  that 
HIVEL2D  accurately  predicted  the  asymmetric  flow  patterns  downstream  of 
the  expansion. 


Supercritical  Confluence 

The  numerical  model  was  further  tested  by  comparing  computed  results  of 
the  Puerto  Nuevo/Guaracanal  Channel  confluence  with  those  observed  in  the 
laboratory  study.  Geometric  details  of  the  supercritical  confluence  are  pro¬ 
vided  in  Figure  12.  Table  3  lists  the  input  parameters,  and  Figure  13  shows 
the  fmite  element  grid,  which  is  compo^  of  491  nodes  and  390  elements. 
Boundary  conditions  for  these  tests  were  supercritical  inflow  with  3.6  cfs 
(0. 10  cms)  in  the  main  channel  and  1 . 1  cfs  (0.03  cms)  in  the  tributary  channel 
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Tabl«3 

Input  P«rain«tar»  for  tho  Suporcritic^  Confkwnco 

a 

1.0 

0.1,06 

n 

0.009 

0.006,0.006  ft'/MC  (4.66  x  lO"  mVsac) 

9 
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and  supercritka]  outflow.  Figures  14  and  IS  show  that  the  HIVEL2D 
captured  the  overall  features  of  the  diamond-shaped  standing  wave  pattern 
resulting  from  the  confluence  geometry.  The  water  surfKc  contours  shown  in 
Figure  IS  illustr^  the  model’s  ability  to  simulate  si4)crelevation  in  the  water 
surface  within  the  tributary  channel  b^.  Water-surface  profiles  along  the 
flume  walls  obtained  using  HIVEL2D  and-those  observed  in  the  laboratory  are 
presented  in  Figure  16.  In  most  siq>ercritical  channels  the  water  surface  oscil¬ 
lates  in  time  even  under  steady  boundary  conditions.  The  recorded  laboratory 
results  are  the  maximum  water-surface  elevations  observed,  whereas  the 
numerical  model  represents  average  water-surface  elevations.  HIVEL2D  ade¬ 
quately  simulated  the  initial  shock  wave  crest,  but  the  location  of  each  subse¬ 
quent  wave  crest  was  increasingly  in  error.  This  difference  is  due  to  the  shal¬ 
low-water  assumption  used  in  the  2-D  numerical  model.  The  shallow-water 
assumption  results  in  all  waves  traveling  with  the  celerity  of  a  long  wave, 
whereas  three-dimensional  flow  is  actually  composed  of  many  wave  speeds, 
the  maximum  of  which  is  the  long-wave  celerity.  The  larger  wave  speed 
means  that  the  standing  wave  angles  will  be  greater  than  the  three-dimensional 
waves.  Generally,  this  shallow-water  equation  limitation  should  be  of  little 
consequence  since  channel  wall  heights  are  set  to  contain  the  maximum  water- 
surface  elevation  plus  freeboard. 


Bridge  Pier 


A  hypothetical  bridge  pier  was  tested  to  demonstrate  HIVEL2D’s  ability  to 
simulate  a  Class  B  bridge  in  which  the  piers  choke  the  flow  such  that  the 
approaching  flow  is  subcritical  even  thtmgh  the  channel  slope  is  hydraulically 
steep.  The  channel  was  100  ft  (30.5  m)  wide  by  1,000  (305  m)  ft  long.  Two 
piers  were  used  to  ctmke  the  flow.  Each  pier  was  20  ft  (6.1  m)  wide  by 
100  ft  (30.5  m)  long  with  a  square  nose  and  tail.  The  finite  element  grid  for 
the  bridge  pier  problem  was  composed  of  993  nodes  and  900  elements 
(Figure  17). 

Water-surface  meslMS  calculated  using  HIVEL2D  are  shown  in  Figures  18 
and  19.  Table  4  lists  input  parameters  for  this  simulation.  The  flow 
approaching  the  piers  was  subcritical  and  accelerated  around  the  piers  to 
supercritical.  The  inqxjsed  tailwater  resulted  in  a  hydraulic  junq)  immediately 
downstream  of  the  piers.  The  hypothetical  problem  demonstrates  HIVEL2D’s 
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T«bl«4 

Input  Parmnotara  for  tha  Brkloa  Plar  Problam 

ConMofi 
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a 

1.0 

AA 

o.z.o.s 

n 

0.020 

v.v. 

10.0,10.0  re/sac  (0.93  m'/aac) 
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ability  to  model  subcritical  flow,  the  transitkm  to  st4)ercritical  flow,  and  also 
the  transition  to  subcritical  flow  via  the  hydraulic  jump. 


Dam  Break  Problem 

Unsteady  flow  sinnilations  were  conducted  to  test  HIVEL2D*s  discretiza¬ 
tion  to  the  temporal  derivatives  of  the  governing  equations.  Results  were 
compared  to  hydraulic  flume  results  repotted  in  Bell,  Elliot,  and  Chaudhry 
(19^).  A  plan  view  of  the  flume  facility  is  shown  in  Figure  20.  The  flume, 
constructed  of  Plexiglas,  simulated  a  dam  break  through  a  horseshoe  bend. 
Because  this  was  a  2-D  problem  and  model  results  were  being  compared  to 
hydraulic  flume  results,  the  limitations  of  the  shallow-water  equations  them¬ 
selves  needed  to  be  considered.  Initially,  the  reservoir  had  an  elevation  of 
0. 1898  m  relative  to  the  channel  bed;  t^  chaiuiel  itself  was  at  a  dq>th  (and 
elevation)  of  0.0762  m.  The  velocity  was  zero  and  then  the  dam  was 
removed.  The  surge  location  and  height  were  recorded  at  several  stations  aixl 
compared  to  the  model  at  three  of  these,  at  stations  4,  6,  and  8.  Station  4  was 
6.00  m  from  the  dam  along  the  channel  center  line  in  the  center  of  the  bend. 
Station  6  was  7.62  m  from  the  dam  near  the  conclusion  of  the  bend,  and 
Station  8  was  9.97  m  from  the  dam  in  a  straight  reach.  The  nnodel-specified 
parameters  are  shown  in  Table  5. 


TablaS 

Test  Conditions  For  The  Dam  Break  Case 

1  Condhion 

VahM 

a 

1 .0,  1 .6 

AA 

0.25,  0.50 

n 

0.009 

v.v. 

0.001,  0.01  m'/sac 

9 

9.802  m/sac’ 

Or 

0.05  sac 
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The  munertcal  grid  (Figure  21)  ccMMiia  698  elements  aynd  81 1  nodes.  This 
grid  was  reached  by  increasing  the  resolutioo  until  the  results  no  longer 
changed.  The  most  critical  reach  was  in  the  region  of  the  contraction  near  the 
dam  breach.  The  basic  element  length  in  the  channel  was  0. 1  m  with  five 
elements  across  the  channel  width.  For  the  smooth  channel  case.  Bell,  Elliot, 
and  Chaudhry  (1992)  used  a  one-dimensional  calculation  to  estimate  the 
Manning’s  n  to  be  0.016,  but  experience  at  the  U.S.  Army  Engineer  Water¬ 
ways  Experiment  Station  suggests  that  this  value  should  actually  be  0.009 
(Brater  and  King  1976). 

The  test  results  for  stations  4,  6,  and  8  are  shown  in  Figures  22-24.  Here 
the  time-history  of  the  water  elevation  is  shown  for  the  inskle  and  outside  of 
the  chaiuiel  for  both  the  numerical  model  (at  a  of  1.0  and  l.S)  and  the  flume. 
The  inside  wall  is  designated  by  squares  aiKi  the  outside  by  diamonds.  Of 
particular  inqwrtance  is  the  arrival  time  of  the  shock  front.  At  station  4  the 
numerical  pr^iction  of  arrival  time  using  a  of  1.0  was  about  3.4  sec,  which 
a{^)ears  to  be  about  0.05  sec  sooner  than  for  the  flume.  This  is  roughly  1- 
2  percent  fast.  For  or  of  1 .5  the  time  of  arrival  was  3.55  sec,  which  is  about 
0. 1  sec  late  (3  percent).  At  station  6  both  flume  and  numerical  mode)  arrival 
times  for  a  of  1.0  were  about  4.3  sec,  and  for  station  8  the  inimerica)  model 
was  5.6  sec  and  the  flume  was  5.65  to  5.8  sec.  With  a  set  at  1.5,  the  time  of 
arrival  was  late  by  about  0.2  and  0. 15  sec  at  stations  6  and  8,  respectively. 
The  flume  at  stations  6  and  8  had  a  earlier  arrival  time  for  the  outer  wave 
compared  to  the  inner  wave.  The  numerical  model  did  not  show  this.  In 
comparing  the  water  elevations  between  the  flume  and  the  numerical  model,  it 
is  apparent  that  the  flume  results  show  a  more  rapid  rise.  The  numerical 
model  is  smeared  somewhat  in  time,  likely  as  a  result  of  the  flrst-order 
temporal  derivative  calculation  of  a  of  1.0.  The  numerical  model  with  a  set 
at  1 .5  shows  an  overshoot.  This  is  likely  a  numerical  artifact  and  not  based 
upon  physics  even  though  this  looks  mu^  like  the  flume  results.  The  surge 
elevations  predicted  by  the  numerical  model  are  fairly  close  if  one  notices  that 
the  initial  elevation  of  the  flume  data  was  supposed  to  be  0.0762  m  and  it 
appears  to  be  recorded  as  much  as  0.015  m  higher  at  some  gauges.  Since  the 
velocity  was  initially  zero,  then  all  of  these  readings  should  have  been 
0.0762  m  and  all  should  be  adjusted  to  match  this  initial  elevation. 

With  this  in  mind,  stations  4  and  8  match  fairly  closely  between  flume  and 
numerical  model.  Station  4  in  the  flume  wiHild  still  have  a  greater  difference 
between  outer  and  inner  wave  than  that  predicted  by  the  model.  The  differ¬ 
ence  might  be  a  manifestation  of  a  threeKlimension^  effect  that  the  model 
cannot  mimic.  The  overall  timing  and  height  conqurisons  are  good. 

Figure  25  shows  the  spatial  profile  of  the  outer  wall  water-surface 
elevation  of  the  numerical  model  versus  distance  downstream  from  center  line 
of  the  dam.  The  two  conditions  are  for  a  of  1.0  and  1.5,  i.e.,  first-  and 
second-order  temporal  derivative. 

The  nodes  are  delineated  by  the  symbols  along  the  lines.  The  overshoot  of 
the  second-order  scheme  and  the  danqiing  of  the  first-order  are  obvious. 
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Again,  it  is  probable  that  the  overshoot  is  a  numerical  artifact  even  though  this 
is  much  like  what  the  flume  would  show. 

Roll  Waves 

Roll  waves  are  a  (dienomenon  of  a  high  Froude  number  environmott.  If 
the  Froude  number  is  greater  than  2.  it  can  be  shown  that  disturbances  can 
grow  into  shocks,  which  progress  downstream  faster  than  the  flow.  This 
results  in  a  characteristic  sawtooth-shaped  water-surface  pattern  in  which  the 
steep  face  is  downstream.  The  roll  wave  tests  used  a  mc^ified  version  of 
HIVEL2D  in  which  a  sinusoidal  perturbation  was  input  at  the  upstream 
boundary  of  a  long  straight  charnel.  The  water  depth  results,  shown  in 
Figure  26,  demonstrate  the  gradual  steepening  of  tte  downstream  wave  face 
downstream  from  inflow  point.  While  this  is  a  qualitative  comparison,  the 
model  does  demonstrate  die  capacity  to  rqiroduce  this  characteristic  event. 


12 


Chaptar  3  Test  CatM 


4  Discussion  and 
Conciusions 


This  series  of  tests  demonstrates  the  abUity  of  the  model  HIVEL2D  to 
supply  engineering  decision  makers  with  a  tool  to  evaluate  hydraulic  results  of 
structural  modificatioos  in  supercritical  channels. 

The  code  itself  is  relatively  flexible  with  the  one  major  limitation  that  the 
boundary  conditions  are  constant  over  a  simulation.  The  most  significant 
limitations  are  therefore  imposed  by  the  equations  modeled.  These  are  the 
shallow-water  equations  employing  the  mild  geometric  slope  assumption  and 
hydrostatic  pressure  distribution.  The  obvious  reailt  is  that  one  would  not 
want  to  use  the  model  to  evaluate  features  with  sleep  slopes.  However,  the 
hydrostatic  assumption  is  more  subtle.  The  result  of  this  assumption  is  that 
shorter  wavelengths  wfll  tend  to  propagate  too  quickly  in  the  model.  This  was 
particularly  noted  in  the  supercritical  confluence  tests.  While  the  water-surface 
predictkMis  were  good  over  the  first  reflection  or  so,  they  became  progressively 
worse  downstreanL  Another  consequence  of  the  hydrostatic  assumption  is  that 
energy  is  dissipated  too  quickly.  The  result  is  that  energy  that  should  be  cap¬ 
tured  in  vertical  motion  is  n^ected  and  undular  jumps  are  instead  modeled  as 
strong  jumps.  Shallow-water  modds  in  which  vertical  motions  are  assumed 
eligible  cannot  predict  the  length  of  a  hydraulic  jump. 
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Figural.  Num«ricai  model  oompuMbonal  mash  for  tha  channel  contraction 
problem 
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Laboratory  flume  (modMed  from  Ippen  wtd  Dawson  1951) 


0.5  0  05  t  o  t  S  2.0  2.5  3.0  15  4.0  4.5  10  15  6.0  6  5  7.0 

OtSTANCt  FROM  COMTAACTXJN.  FT 


b.  Numerical  modal  (HIVEL20) 


Figure  2.  Depth  contours  of  supercritical  flow  in  a  contraction 


Figure  3.  Water-surface  mesh  of  supercritical  9am  in  a  corUraction 


Figur*  4.  Margwla  Channai  fluma  datatfs  Figura  5.  Numarical  modal  oomputationai  mash  (or  tha  Margarita 


Figure  6.  Water- surface  profiles  of  Margarita  Channel  flume 


Figure  7.  Lateral  dbtribiiliortt  of  kmgitiKinai  velocity  in  Margarita  Channel  at 
station  1 


Figure  8.  Lateral  distribulions  of  longNudirud  velocity  in  Figure  9.  Lateral  distributio 

Margarita  Channel  at  station  2  velocity  in  Margai 

station  3 


Figure  1 1 .  HIVEL2D  results  of  Margarita  Channel.  Velocity  contours  and 
vectors  (darker  contours  indicate  lower  velocities) 


1.714-FT  VWDTH — ■ 
0.000  FT 
MODEL  LIMIT 


Figure  12.  Details  of  supercritical  confluence 


CHANNEL 


Figure  13.  Numerical  model  computational  mesh  for  the  supercritical 
confluence  problem 


a.  Right  side 


b.  Leftside 


Figure  16.  Water-surface  profiles  for  the  supercritical  confluence  problem 


Figure  17.  Numerical  model  computational  mesh  for  the  bridge  pier  problem 


Figure  19.  Side  view  of  the  water-surface  mesh  of  the  bridge  pier  problem 
(flow  is  from  right  to  left) 


Figure  20.  Details  of  the  dam  break  test  flume  (modified  from  Bell,  Elliot,  and 
Chaudhry  1992) 


Figure  22.  Flume  and  numerical  model  depth  histories  for  station  4 
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a.  Fhime 


c.  Numerical  model,  a,  ■  1.5 


Figure  24.  Fhjme  arKi  ruimerical  model  depth  Nstories  for  station  8 


Figure  26.  Depth  profiles  for  the  roll  wave  probienn 


Appendix  A 
Governing  Equations 


Vertical  integration  of  the  three-dimensional  equations  of  mass  and  momen¬ 
tum  conservation  for  incompressible  flow  with  the  assumption  that  vertical 
velocities  and  accelerations  are  n^igible  compared  to  horizontal  motions  and 
the  acceleration  of  gravity  results  in  the  governing  equations  commonly 
referred  to  as  the  shallow-water  equatkms.  The  dcpeadent  variables  of  the 
two-dimensional  fluid  motion  are  defined  by  the  flow  dq>th  h,  the  x-direction 
component  of  unit  discharge  p,  and  the  y-direction  component  of  unit  dis¬ 
charge  q.  These  variables  are  functions  of  the  indqrendent  variables  x  and  y, 
the  two  ^ace  directions,  and  time  t.  Neglecting  fiee-surface  stresses  and  the 
effects  of  Coriolis  force  as  these  are  not  considered  important  in  high-vdocity 
chatmels,  the  shailow-water  equations  in  conservative  form  are  given  as 
(Abbott  1979;  Praagman  1979):^ 


♦  ii  -  0  (Al) 

dt  dx  dy 


for  the  conservation  of  mass.  Conservatktn  of  momentum  in  the  x-direction 
and  y-direction  are  given  res^rectively  as: 


+  -i 

dt  h 


4- 


2 


(A2) 


and 


*  Reterenoes  died  in  this  Appendix  are  listed  in  Ihe  References  at  the  end  of  the  main  text. 
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Al 


(A3) 


h  *  A.\es. 

dt  dx  h 


Ifi! 

3y  h 


n^q^/p^ 


ho„ 


where 


g  =  acceleration  of  gravity 

a„.a^o^,a^  =  Reynolds  stresses  per  unit  mass  where  the  first  subscript 
indicates  the  direction  and  the  second  indicates  the  face  on 
which  the  stress  acts 

z  -  channel  invert  elevation 

n  =  Maiming’s  roughness  coefficient 

Q  =  dimensional  constant  1  for  SI  units  and  2.208  for 
non-SI  units) 

The  governing  equations  are  given  in  vector  form  as: 


dt  da 


dF 

+  _I  +  H  =  0 
dy 


(A4) 


where 


Q  = 


h 

P 

q 


(A5) 


(A6) 
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where 


ehii  - »» gyp  q 

CoV/3 


I 


pss  uh,u  being  the  depth-averaged  x-direction  component  of  velocity 

q  s  vh,v  being  the  d^th-averaged  y-direction  component  of  velocity 

The  Reynolds  stresses  are  determined  using  the  Boussinesq  approach  of 
gradient-diffusion: 


(A9) 


Where  v,  is  the  viscosity  (sum  of  turbulent  and  molecular  viscosity,  commonly 
refened  to  as  eddy  viscosity),  which  varies  ^atially  and  is  solved  empirically 
as  a  function  of  local  flow  variables  (Rodi  1980;  Chapman  and  Kuo  1985): 


V 


t 


*1/6 


(AlO) 


where  is  a  coefHcient  that  varies  betwera  0.1  and  1.0. 
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This  system  of  equations  constitutes  a  hyperbolic  initial  boundary  value 
problem.  Appropriate  boundary  conditions  are  determined  using  the  approach 
of  Daubert  and  Graffe  as  discussed  in  Drolet  and  Gray  (1988)  and  Vert^m, 
Stelling,  and  Officier  (1982).  Daubert  and  Graffe  use  the  method  of  charac¬ 
teristics  to  determine  the  required  boundary  conditions.  The  number  of  bound¬ 
ary  conditions  is  equal  to  the  number  of  characteristic  half-planes  that  originate 
exterior  to  the  domain  and  enter  it  If  the  inflow  boundary  is  supercritical, 
then  all  information  from  outside  the  domain  is  carried  through  this  boundary. 
Therefore,  p  and  q  (or  u  and  v)  and  the  dqjth  h  must  be  specified.  If  the 
inflow  boundary  is  subcritical,  then  the  depth  is  influenced  from  the  flow 
inside  the  domain  (downstream  control)  and  therefore  only  p  and  q  (or  u  and 
v)  are  specifled.  Outflow  boundary  conditions  required  are  determined  by 
analysis  of  information  transported  through  this  boundary.  If  the  outflow 
boundary  is  supercritical,  then  all  information  is  determined  within  the  domain 
and  no  boundary  conditions  are  specified.  However,  if  the  outflow  boundary 
is  subcritical,  then  the  depth  of  flow  at  the  boundary  (tailwater)  must  be 
specifled.  The  no-flux  boundary  condition  is  appropriate  at  the  sidewall 
boundaries  and  is  discussed  in  detail  in  Appendix  B. 
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Appendix  B 

Finite  Element  Formulation 


A  variational  formulation  of  the  governing  equations  involves  finding  a 
solution  of  the  dependent  variables  Q  using  the  test  function  over  the 
domain  0.  The  variational  formulation  of  the  shallow-water  equations  in 
integral  form  is: 


dt 


dF.  ap. 


dx 


dy 


^  H 


dO  =  0 


(Bl) 


where  t  is  time  and  (2^  F,  and  iY  are  defined  in  Equations  AS-A8. 

The  finite  element  approach  taken  is  a  P^rov-Galerldn  formulation  that 
incorporates  a  combination  of  the  Galerkin  test  function  and  a  non-Galerldn 
component  to  control  oscillations  due  to  convection.  The  finite  element  form 
of  the  governing  equations  is: 


E 


■ 

f 

.  fl 

dO, 

L  4 

dt  dx  dy 

€ 

=  0,  for  each  i 

(B2) 


where 

e  =  subscript  indicating  a  particular  elenwnt 
i  —  subscript  indicating  a  particular  test  function 
~  =  discrete  value  of  the  quantity 

The  geometry  and  flow  variables  ate  represented  using  the  Lagrange  basis 
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Bl 


(B3) 


1?  =  E 

j 


where  y  is  the  nodal  location.  Bilinear  triangular  and  quadrilateral  elements 
are  us^  with  nodes  at  the  element  comers.  Figure  B1  shows  the  two  bilinear 
elements  used  in  terms  of  local  coordinates  (  and  i}. 


Figure  B1.  Local  bilinear  elements 


The  test  function  used  (to  be  elaborated  in  the  mxt  section)  is: 


(B4) 


where 

0,  =  Galerkin  part  of  the  test  function 
/  =  identity  matrix 

=  non-Galerkin  part  of  the  test  function 

To  facilitate  the  specification  of  boundary  conditions,  the  weak  form  of  the 
equations  is  develq)^  using  an  integration  by  parts  procedure.  Integration  by 
parts  of  the  terms 


B2 
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yields  the  weak  form  of  the  equations.  The  is  omitted  for  clarity  and  the 
variables  are  understood  to  be  discrete  values.  The  weak  form  is  given  as: 


j  [  ‘  dt  3x  ^  dy  ^  ‘  dx  ‘  dy 

dO,  +  f<l>.  >  Fyi,)  dT, 


(B6) 


where  (n,,  n^)  =  the  unit  vector  outward  normal  to  the  boundary  F,,  and 


A 

B 


dF, 

dQ 

dQ 


(B7) 


Natural  boundary  conditions  are  tqTplied  to  the  sidewall  boundaries  through  the 
weak  statement.  The  sidewall  boundaries  are  ”no  flux”  boundaries.  That  is, 
there  is  no  net  flux  of  mass  nor  momentum  through  these  boundaries.  This 
boundary  condition  is  enforced  in  an  avenge  sense  through  the  weak  state¬ 
ment.  Setting  the  mass  flux  through  the  sidewall  boundary  to  zero: 

f  (pn,  +  gn^  dT  -  0  (B8) 

r 


where 

p  =  x-direction  conqMnent  of  unit  discharge 

q  =  y-direction  conqwnent  of  unit  discharge 

There  is  no  net  momentum  flux  through  the  boundaries.  Therefore,  the 
x-direction  momentum  through  the  boundary  is  set  to  zero: 

Appendix  B  Finita  Element  Formulation 


(B9) 


j  [(»fpK  +  =  0 

r 


and  the  y-direction  momentum  through  the  boundary  is  set  to  zero: 


j  pP)n^  *  (yq)n^  dT  =  0  (BIO) 

r 


where 

u  =  p/h  =  Hot  depth-averaged  x-direction  component  of  velocity 
V  =  q/h  =  ibt  dq>th-averaged  y-direction  conq)onent  of  velocity 
A  =  the  depth  of  flow 

Sidewall  drag  is  treated  as  a  partial  slip  condition.  That  is  the  boundary 
stress  terms  in  the  governing  equations,  integrated  along  the  sidewall,  are 
specified  via  the  Manning  relation: 

-  I  *  ho^y)<^  =  I 

and 

-  I  {hOy/i,  +  ha^y)  rfr  =  I  gq!^:!E^L  dI^12) 

where 

a^a^ay^o„  =  Reynolds  stresses  per  unit  mass  where  the  first  subscript 
indicates  the  direction  and  the  second  indicates  the  face 
on  which  the  stress  acts 

g  —  acceleration  of  gravity 

Cg  =  dimensional  constant  (Co=  1  for  SI  units  and  2.208  for 
non-SI  imits) 
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Petrov-Galerkin  Test  Function 


For  the  shallow-water  equations  in  conservative  form  (Equation  B2),  the 
Petrov-Galerldn  test  function  p,  is  defined  as  (Berger  1993)': 


•Pi  ^  & 


A  i  A 

Ax-^i4  +  Ay-:-LB 
Bx  'By 


(B13) 


where  is  a  dimensionless  number  b^een  0  and  0.5,  and  ^  is  the  linear 
basis  function.  In  the  manner  of  Katopodes  (1986)  the  grid  intervals  are 
chosen  as: 


Aj:=2 


— 

■"N 

a* 

2 

\^V 

Ll«J 

+ 

dljj 

(B14) 


and 


Ay  =2 


(B14) 


where  {  and  i;  are  the  local  coordinates  defiiKd  from  -1  to  1  (Figure  Bl). 


To  find  A  consider  the  following: 

A  m 

(B16) 

BQ 

P-^  AP  =  A 

(BIT) 

where  A  =  A  is  the  matrix  of  eigenvalues  of  A,  and  P  and  are  made  up 
of  the  right  and  left  eigenvectors. 


'  References  cited  in  this  tppendix  are  listed  in  die  References  at  the  end  of  the  main  text. 
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Am  P-^  A  P 


(B18) 


where 


+  V^f 


A 

A  = 


0 


0 


(B19) 


and 


X,  =  u  +  c 

(B20) 

Xj  =  It  -  c 

(B21) 

Xj  =  M 

(B22) 

c  =  (gh)''^ 

(B23) 

A  similar  operation  may  be  performed  to  define  A 

This  particular  test  function  is  weighted  upstream  along  characteristics 
similar  to  a  concept  like  that  developed  in  the  finite  difference  method  of 
Courant,  Isaacson,  and  Rees  (1952)  for  one-sided  differences.  These  ideas 
were  expanded  to  more  general  problems  by  Moretti  (1979)  and  Gabutti 
(1983)  as  split-coefficient  matrix  methods  and  by  the  generalized  flux  vector 
splitting  proposed  by  Steger  and  Warming  (1981).  In  the  finite  element 
community,  instead  of  one-sided  differences  the  test  function  is  weighted 
upstream.  This  particular  method  in  one  dimension  (1-D)  is  equivalent  to  the 
SUPG  (streamline  upwind  Petrov-Qalerkin)  scheme  of  Hughes  and  Brooks 
(1982)  and  similar  to  the  form  proposed  by  Dendy  (1974).  Examples  of  this 
approach  in  the  open  channel  environment  usii^  Ae  generalized  shallow-water 
equations  are  presented  for  1-D  in  Berger  and  Winant  (1991)  and  for  2-D  in 
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Berger  (1992).  A  1-D  St.  Venant  application  is  given  by  Hicks  and  Steffler 
(1992). 


Shock  Capturing 

Berger  (1993)  shows  that  the  Petrov-Galerkin  scheme  is  not  only  a  good 
scheme  for  advection-dominated  flow,  but  is  also  a  good  scheme  for  shock 
capturing  because  the  scheme  dissipates  energy  at  the  short  wavelengths. 

When  a  shock  is  encountered,  the  weak  solution  of  the  shallow-water 
equations  must  lose  mechanical  energy.  Some  of  this  energy  loss  is  analogous 
to  a  physical  hydraulic  system  losing  energy  to  heat,  particle  rotation,  etc;  but 
much  of  it  is,  in  fact,  sin:^)ly  the  energy  being  transferred  into  vertical  motion. 
And  since  vertical  motion  is  not  included  in  the  shallow-water  equations,  it  is 
lost.  This  apparent  energy  loss  can  be  advantageous. 

To  apply  a  high  value  of  0,  say  O.S,  only  in  regions  in  which  it  is  needed, 
since  a  lower  value  is  more  precise,  construct  a  trigger  mechanism  that  can 
detect  shocks  and  increase  0  automatically.  The  method  employed  detects 
energy  variation  for  each  element  and  flags  those  elements  that  have  a  high 
variation  as  needing  a  larger  value  of  0  for  shock  cq)turing.  Note  that  this 
would  work  even  in  a  Galerkin  scheme  since  this  trigger  is  concerned  with 
energy  variation  on  an  element  basis  and  the  Galerkin  method  would  enforce 
energy  conservation  over  a  test  function  (which  includes  several  elements). 

The  shock  capturing  is  implemented  when  Equation  B24  is  true 

Tr,  >  7  (B24) 


where  7  is  a  specified  constant  and 


Ts,  = 


£D,- 


(B25) 


where  £Z>„  the  element  energy  deviation,  is  calculated  by 


EDi  = 


(E  -  E^  dQ 


la 


(B26) 
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where 


0^  =  element  < 

E  =  mechanical  energy 
n,  =  area  of  element  i 

and  £ the  average  energy  of  elonent  i,  is  calculated  by 


Eda 


n. 


(B27) 


and  _ 

£  =  the  average  element  energy  over  the  entire  grid 
5  =  the  standard  deviation  of  all  £D, 

Through  trial  a  value  of  7  of  1.0  was  chosen. 


Temporal  Derivatives 

A  finite  difference  expression  is  used  for  the  temporal  derivatives.  The 
general  expression  for  the  temporal  derivative  of  a  variable,  Qj,  is: 


9Qj 

m*\ 

»  a 

fiT  '  -  (T 

+  (1  -  a) 

-cT' 

dt 

where 

a  =  temporal  difference  coefficiait 
j  =  nodal  location 


m  =  time-stq) 


An  a  equal  to  1  results  in  a  first-order  backward  difference  approximation  and 
an  a  equal  to  1.5  results  in  a  second-order  backward  difference  a{^roximation 
of  the  temporal  derivative. 
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Solution  of  the  Nonlinear  Equations 


The  system  of  nonlinear  equations  is  solved  using  the  Newton-Raphson 
iterative  method  (Carnahan,  Luther,  and  Wilkes  1969).  Let  be  a  vector  of 
the  nonlinear  equations  computed  using  a  particular  test  function  ^id  using 
an  assumed  value  of  R,  is  the  residual  error  for  a  particular  test  function 
i.  Subsequently,  Ri  is  forced  toward  zero  as: 


^  =  -R*  (B29) 

ac* 


where  k  is  the  iteration  number,  j  is  the  node  location,  and  the  derivatives 
composing  the  Jacobian  are  determined  analytically.  This  system  of  equations 
is  solved  for  Aq/  and  then  ai:  improved  estimate  for  is  obtained  from: 


a;'  A,* 


(B30) 


This  procedure  is  continued  until  convergence  to  an  accej^le  residual  error 
is  obtained. 

Equation  B29  represents  a  system  of  linear  algd}raic  equations  that  must  be 
solved  for  each  iteration  and  e^  time-step.  A  ■Profile"  solver  is  imple¬ 
mented  to  achieve  efficient  coefficient  matrix  storage.  This  method  stores  the 
upper  triangular  portion  of  the  coefficient  matrix  by  columns  and  the  lower  by 
rows.  Any  zeros  outside  the  profile  arc  not  stored  or  involved  in  the  com¬ 
putations.  The  necessary  arrays  are  then  a  vector  composed  of  the  columns  of 
the  upper  portion  and  a  pointer  vector  to  locate  the  diagonal  entries. 

Triangular  deconqrasition  of  the  coefficient  matrix  is  used  in  a  direct  solution. 
The  program  to  construct  the  triangular  decomposition  of  the  coefficient  ma¬ 
trix  uses  a  compact  Crout  variation  of  Gauss  Elimination. 
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